Superconductivity with the Meron-Cluster Algorithm 

J.C. Osborn a * 

department of Physics, Box 90305, Duke University, Durham, NC 27708, U.S.A. 

The meron-cluster algorithm was previously used to extensively study the physics associated with the sponta- 
neous breaking of a discrete symmetry. We recently discovered that a larger class of models with spontaneous 
breaking of continuous symmetries can also be simulated using the meron-cluster algorithm. Here we study one of 
these new models which belongs to the attractive Hubbard model family. In particular we study the spontaneous 
breaking of the (7(1) fermion number symmetry in two dimensions and find clear evidence for a Kosterlitz-Thouless 
transition to a superconducting phase. 



1. Introduction 

For over a decade there has been a lot of 
interest in performing numerical simulations of 
strongly correlated fermionic systems. A recently 
discovered method for solving the sign problem 
inherent in such simulations is the meron-cluster 
algorithm A new reference configuration al- 
lows us to explore a larger class of models with 
this method. The model we study here is a vari- 
ant of the attractive Hubbard model. 

The attractive Hubbard model has recieved a 
lot of attention as a toy model for superconduc- 
tivity. It is a very simplistic model since the at- 
traction between fermions is explicitly included, 
instead of arising from some more subtle mecha- 
nism. In that respect it is perhaps more closely 
related to the case of color superconductivity in 
QCD. It is still useful as a model for high-T c 
superconductors since at intermediate attraction 
strength they both exhibit pairing at tempera- 
tures above the superconducting transition || . 

Despite the Hubbard model's simplicity, nu- 
merical attempts to verify the properties of a 
Kosterlitz-Thouless (KT) transition are inconsis- 
tent J3|. The main reason seems to be the fairly 
small lattice size restriction imposed by conven- 
tional fermionic Monte Carlo methods. Since 
the meron-cluster algorithm scales very efficiently 
with lattice volume, we are able to go to large 
system sizes needed to accurately test finite size 
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scaling formulas. Another advantage comes from 
performing the simulation directly in the fermion 
occupation basis, thus allowing very simple access 
to a wide range of observables. 

2. The Model 

We study a model of fermions with spin on a 
lattice in two spatial dimensions. The Hamilton 
operator can be written as 



H — ^ h xy + ^ h x , 

<x,y> x 

with a nearest neighbor interaction 

hxy = ^ t ( c x,s c y,s + c y,s c x.s) 
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Figure 1. Pair correlation versus temperature for 
C/ = 8 and A = 1. 
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Figure 2. Pair correlation versus temperature for 
U = and A = 1. 



along with the on-site spin 

&x — 2 c ^' s <Ts ' s ' C;E:S ' ^ 
and pseudo-spin operators 

Jl = ^ [4+4-1] . (6) 

The term containing U serves to bind the 
fermions with opposite spin together on a lattice 
site. In the limit as U — > oo the model becomes 
bosonic and can be mapped onto the anisotropic 
quantum Heisenberg model. As U decreases the 
model becomes increasingly fermionic, as in the 
attractive Hubbard model. However, at U = 
we still have a strongly interacting system due to 
the additional terms present. 

In spite of these additional factors, we still ex- 
pect this model to behave similarly to the attrac- 
tive Hubbard model since they both possess the 
same symmetries. The important symmetries are 
the SU(2) S spin group generated by the stan- 
dard fermionic spin operator S = ^2 X S X , and 
the SU(2) C charge symmetry generated by the 
pseudo-spin operator J = J x . The SU(2) S 
group is always a symmetry of the Hamiltonian, 
but the SU(2) C group is broken down to the 
U(\) c particle number group if either /i or A 
are nonzero. The remaining U{l) c symmetry can 
then undergo a Kosterlitz-Thouless transition in 
two dimensions. 



Here we will only consider the case where the 
SU(2) C is broken by A ^ at half filling (fj, = 0) 
as a test of the model and our simulations. We 
will study the more physically relevant case of 
fj, 7^ elsewhere ||]. 

3. Simulation Method 

The meron-cluster algorithm uses clusters in 
the standard path integral formulation of a quan- 
tum model to solve the sign problem arising from 
fermion permutations. It has been presented in 
detail for a variety of models . Here we will 

only outline the extension to this model. 

We chose this particular model because of its 
simple cluster rules. The clusters are formed iden- 
tically for the spin up and spin down fermions. 
We can thus group all the clusters into pairs 
that are identical except for having different spin. 
Flipping both clusters to be in the same orienta- 
tion then gives an identical fermion sign for each 
spin sector, thus canceling each other out. No 
matter how the clusters are shaped, we can al- 
ways flip them to be in a "reference configura- 
tion" with both spin sectors identical, which will 
always have a positive contribution to the parti- 
tion function. 

The on-site attraction term proportional to U 
is implemented with an extra on-site plaquette 
similar to the mass term in the previously stud- 
ied staggered fermion model [||. Here, however, 
the extra link can either do nothing, or it can 
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Figure 3. Winding number versus temperature 
for Z7 = 8 and A = 1. 
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Figure 4. Winding number versus temperature 
for U = and A = 1. 



freeze two clusters of different spin together in the 
same orientation. If frozen together the clusters 
can still be flipped, but both spin sectors must 
be flipped together. Note that this frozen state is 
consistent with the reference configuration which 
allows us to include it in the meron-cluster algo- 
rithm. This also limits us to studying an on-site 
attraction in this model since a repulsion would 
need to freeze the clusters in opposite orientations 
which is contrary to the reference configuration. 

4. Numerical Results 

We simulated the model presented above using 
a variant of the implementation used in HQ . We 
have been able to simulate systems with a spa- 
tial volume (V) of 128 2 , but we will only present 
results of up to 64 2 here. Each simulation pro- 
duced between 10 5 and 10 6 configurations after 
performing at least 10 thermalization sweeps. In 
all cases we measured the autocorrelation time to 
be less than two configurations. 

The transition to a state with long range order 
can be seen in the S-wave pair correlation 
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In the thermodynamic limit (L — > oo), the pair 



correlation should diverge near a KT transition 
as oc £(^/ 4 with 
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Figures |l| and |^ show the pair correlation for 
U = 8 and U — respectively. In both cases we 
can clearly see evidence of long range correlations 
forming. Here we do not see much of a difference 
between the nearly bosonic (U — 8) and fermionic 
(U = 0) model, only a slight shift in the critical 
temperature. 

One could attempt to extract T c from fitting 
the pair correlation for the largest volume avail- 
able to the above form for P^. This generally 
turns out to be difficult due to the large volumes 
necessary to approximate Poo near T c . 

An easier way of measuring T c comes from the 
helicity modulus HH which can be defined in 
terms of the winding number as 



2 x x 
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where W x (W y ) is the total number of particles 
winding around the boundary in the x (y) direc- 
tion. This is very convenient to work with since 
we know the finite size scaling form to be 
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with A(T C ) — 1. In figures || and |J we show the 
quantity ttT /T for the two different values of U. 




Figure 5. Pairing fraction versus reduced temper- 
ature for (U = 8) and fermionic (U = 0) limits. 

Here we can clearly see the tendency of a uni- 
versal jump between and 2 at T c . Near T c we 
can fit the helicity modulus for several volumes to 
Ji"l| ) to get the two parameters A(T) and L (T). 
We see that the coefficient A(T) moves approxi- 
mately linearly through 1 as T goes through T c . 
By fitting a straight line to A(T) we can deter- 
mine T c . For U = 8 we estimate that T c = 1 .39(2) 
while for U = we get T c = 1.28(2). 

For A = f with U — » oo our model maps 
onto the quantum XY model with an extra factor 
of 1/4 multiplying the temperature. Scaling the 
measured result for the XY model ||, we should 
get T c — 1.3708(2) for large enough U. Our result 
for U — 8 is consistent with this which provides a 
good check of our simulations. We also see that 
U = 8 is indeed close to the bosonic limit. 

We can better understand how the model 
changes as we adjust U by measuring the pair- 
ing fraction which we define as 



(12) 



In figure g we show the pairing fraction for both 
U = 8 and U = versus the reduced temper- 
ature t = (T — T c )/T c using T c as determined 
above. Here we can distinctly see the effects of 
the on-site attraction. For U — 8 the pairing 
fraction gradually increases as T c is approached 
from above and comes within a few percent of its 
value at T c well before the KT transition. For 
U = we see a more dramatic rise as T c is ap- 



proached from above. Here the pairing fraction 
levels off much closer to T c than for U = 8, but 
still just before the transition. 

5. Conclusions 

We have efficiently simulated a fermionic model 
with the same symmetries as the attractive Hub- 
bard model. In this model we observed a KT 
transition at a finite temperature for A = 1 in 
the bosonic (U=8) and fermionic (U=0) limits. 
By measuring the helicity modulus, we were able 
to easily obtain the critical temperature. 

In the future we will include a chemical poten- 
tial in our simulations. This will allow us to move 
away from half filling. We also plan to study this 
model in 3 dimensions to observe the physics of 
massless Goldstone bosons. 
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